Ballistic missile defense has become a significant priority as national intelligence indicates a growing missile threat from rogue nations that might obtain weapons of mass destruction and use ballistic missiles to fire them at U.S. forces abroad, U.S. allies or the continental United States. A desirable engagement strategy against ballistic missiles is to intercept the target as early as possible during the boost phase or early ascent phase when the target is a large object and has not dispersed countermeasures or multiple warheads. Such a strategy minimizes the requirements for warhead and decoy discrimination, and allows for multi-layered defense opportunities. A missile defense system supporting this strategy must include an accurate boost phase target state estimator. Without accurate target state estimates, a fire control system cannot obtain valid intercept solutions for launching the interceptor during boost, and intercepting during the early ascent phase.
Challenges in developing a boost phase tracking algorithm include unpredictable target accelerations while the target is in powered flight, and uncertainty of its burn time. Targets powered by solid rocket motor targets present the additional challenge of irregular thrust acceleration profiles. Due to the significant changes in acceleration during the target's boost phase, filters that assume constant acceleration cannot meet the stringent accuracies required of most fire control systems. Current state-of-the-art template-based filters use position and velocity templates assuming constant accelerations or rocket equation acceleration modeling. Such templates are subject to error attributable to motor burn variations, energy management (lofted/depressed trajectories), ISP variations, and early thrust terminations. FIG. 1 is a plot of target thrust versus time for “hot,” nominal and “cold” rocket motors or engines. The nominal thrust profile as a function of time after launch (TAL) is shown by plot 12. The actual motor may be a variation about this nominal profile. For a hot motor, the motor gain is greater than 1, K>1, with the thrust profile shown by plot 14. For a cold motor, the motor gain K is less than 1, K<1, with the thrust profile shown by plot 16. It will be clear that the accelerations represented by the thrust profiles of FIG. 1 are significantly different.
In addition, uncertain knowledge of the time of launch of the target results in error in estimating the time after launch. The uncertainty in knowledge of the time after launch leads to error in indexing into the template, which in turn tends to introduce more error in the acceleration estimate. For example, if the estimate of time after launch in FIG. 1 is indicated by testimate and the actual time is indicated by tactual, there will be a substantial error in estimating the acceleration even if the proper plot were known. The indexing error and acceleration variations from the nominal template all contribute to erroneous acceleration estimates and ultimately poor estimates of the target velocity and position. The indexing error, delT, is the uncertainty in the target's time after launch (TAL), i.e., the time difference between the initial index into the nominal template and the target's true TAL.
U.S. patent application Ser. No. 11/189,234 filed Jul. 26, 2005 in the name of Luu et al. and entitled Template Updated Boost Algorithm (TUBA) describes a method for estimating the states of a boosting target using nominal templates. The templates include profiles of the target (a) thrust acceleration, (b) altitude, (c) speed, and (d) angle-of-attack as a function of time after launch, all updated using the states of a filter to correct for template indexing error, template boost acceleration variations, and template angle-of-attack variations, thereby producing estimated target states unencumbered by motor burn boost variations and providing for accurate determination of end-of-boost.
According to another aspect of the method for estimating at least the position and velocity state of a boost vehicle target according to the Luu et al. application, the state estimates include the position and velocity, and at least one of (a) variation in engine burn rate, (b) error in time after launch and (c) error in angle of attack. The method includes the step of sensing at least the position of the target to thereby generate target-position information. The method also includes the step of providing a template characterizing nominal values of acceleration of the boost vehicle target due to engine thrust, speed, altitude, and angle-of-attack, all as a function of time after target launch. The target altitude is calculated or determined, possibly from the target-position information, to generate calculated target altitude. A time index is established by entering the template at the calculated target altitude. At least that (or those) nominal values of time after launch, acceleration, speed, and angle of attack corresponding to the calculated target altitude are read from the template. Both the state estimate and the associated covariances of an extended Kalman filter are initialized with the nominal value of time after launch, acceleration, speed, and angle of attack corresponding to the calculated target altitude. The state estimates and covariances of the Kalman filter are propagated in time to generate a time-updated vector of state estimates representing the predicted state of the target, using the template profile information updated by at least one of (a) variation in engine burn rate, (b) error in time after launch and (c) error in angle of attack. The states are updated using the measurements, to thereby produce a vector of measurement updated state estimates and covariances of the boost vehicle target, including position and velocity, and at least one of (a) variation in engine burn rate, (b) error in time after launch and (c) error in angle of attack.
The Template Updated Boost Algorithm (TUBA) is a boost phase filter or processing that estimates variations with respect to nominal templates representing the target's kinematic motion. FIG. 2 is an example of a typical nominal template file of thrust acceleration versus time in seconds (s). In FIG. 2, the thrust acceleration profile of a solid-fuel motor is illustrated by a solid line designated 210, and the profile for a liquid-fuel motor is illustrated by a dash line designated 212. TUBA is capable of estimating rocket motor burn variations, template indexing error (the time difference between the actual and the estimated Time After Launch), and the target's angle of attack. Using nominal templates corrected by these unique filter states, TUBA predicts future thrust acceleration, position, and velocity vectors with improved reliability. This can be crucial for fire control systems, where knowing the target's future position and velocity is needed to accurately predict an intercept point.
Another boost phase filter is described in U.S. Pat. No. 7,181,323, issued Feb. 20, 2007 in the name of Boka et al. and entitled Computerized Method for Generating Low-Bias Estimates of Position of a Vehicle From Sensor Data. This boost filter is the Unified Unbiased Rocket Equation Extended Kalman Algorithm (UUREEKA). It differs from the current TUBA in that UUREEKA models target dynamics using the fundamental rocket equation and is ideal for tracking liquid-fueled targets whose acceleration profiles can be modeled using the rocket equation. UUREEKA is less advantageous for solid fuel rocket motor that exhibit irregular thrust profiles. TUBA is well adapted for tracking solid-fuel targets that have irregular thrust acceleration profiles, which cannot be modeled by the rocket equation. As mentioned, FIG. 2 shows a representative thrust acceleration profile 210 for a solid rocket motor, and a liquid rocket motor profile 212. Since TUBA is template based, it can also be used to track liquid-fuel targets governed by the rocket equation. TUBA complements UUREEKA by providing capability against solid rocket fuel based target threats. Thus, the combination of UUREEKA and TUBA filtering may provide more complete capability against all ballistic missile threats.
TUBA is a template filter which utilizes a table of nominal target data relating target time after launch (TAL) to boost acceleration, speed, altitude, and angle of attack. An example of a nominal target template data table 300 appears as FIG. 3. In FIG. 3, various Time After Launch (TAL) values designated T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, and T11 correlate to Boost Acceleration values AB1, AB2, AB3, AB4, AB5, AB6, AB7, AB8, AB9, AB10, and AB11, respectively, to Speed values V1, V2, V3, V4, V5, V6, V7, V8, V9, V10, and V11, respectively, to Altitude values Alt1, Alt2, Alt3, Alt4, Alt5, Alt6, Alt7, Alt8, Alt9, Alt10, and Alt11, respectively, and to Angle of Attack values AOA1, AOA2, AOA3, AOA4, AOA5, AOA6, AOA7, AOA8, AOA9, AOA10, and AOA11, respectively. Of course, a template may include many more values than the illustrated 11 values, or possibly fewer.
FIG. 4 shows how the nominal target template is used in the TUBA process. In FIG. 4, the current radar target measurement is applied by way of a path 410. Ordinarily, the template initialization is performed using the altitude template. Systems such as Spaced Based Infrared Sensor (SBIRS) also provide velocity information which could be used for template initialization using the speed template.
The entry point into the target template 300 of FIG. 4 is found using the current measurement altitude, Alto, calculated from current target measurement data. A lookup algorithm is used to select the time index, Io, in the table corresponding to Alto. Io refers to the template time index TI indicated by a subscript in the TAL column of FIG. 3. Thus, if the calculated (or otherwise determined) altitude is, say, one thousand feet, the corresponding altitude in the table might be T4, and the corresponding index is I4. Io is then used as the initial time entry into the acceleration and angle of attack templates to obtain the target's nominal boost acceleration, and angle of attack. In the example, the index I4 would point to acceleration AB4 and angle of attack AOA4. TUBA processing as represented by block 412 in FIG. 4 uses the Boost acceleration AB(K, Im) and Angle of Attack AOA(Im) from the template 300 to produce accurate estimates of the target's position P(tm) and velocity states V(tm). In addition, block 412 estimates the error de/A(tm) in angle of attack, error delT in the time index into the template, and a scalar (hot/cold) motor variation factor, K. In conjunction with the typical position and velocity states, additional states delT and K, are used to adjust the initial time index Io to account for initial time index error (possibly attributable to incorrect initial estimate of altitude) and template difference due to motor variations. The TUBA processing produces estimates which are coupled for use by way of paths designated together as 414.
The TUBA algorithm associated with block 412 of FIG. 4 is represented by a logic flow chart or diagram 500 of FIG. 5. In the logic flow of FIG. 5, a measurement is received at time tm as suggested by logic path 510. The algorithm begins with initialization of filter and timing parameters in blocks 512 and 514. Initial estimates of the filter states and state covariances are then propagated forward in time corresponding to the Kalman “time update” step 516. Staging timing parameters are updated in block 524, after the state propagation, and staging events and early thrust termination are determined in block 542. Finally, the Kalman Gains are computed, and the measurement updated state covariances are updated corresponding to the Kalman filter “measurement update” step 526. The process repeats with the measurement updated state estimates and error covariance feeding back to the time update equations.
More specifically, the TUBA process of FIG. 5 begins with arrival over a path 510 of an initial position measurement from a sensor tracking a boosting target. Decision block 512 initially routes the measurement to the TUBA Initialization function or block 514. The TUBA Initialization function 514 initializes the filter and template parameters. FIG. 6 represents the logic 600 of the TUBA Initialization step 514 of FIG. 5. In FIG. 6, the first step in initialization of the filter and template parameters is to calculate target altitude Alto from the current measured position input, as suggested by block 610 of FIG. 6. From block 610, the logic of the initialization 600 flows to a block 612, which represents use of the table lookup to determine target Time After Launch (TAL) template_index, for the given Alto. Block 614 represents use of the template_index produced by block 612 together with nominal staging times. Nominal staging times are not derived from part of the templates which are a function of time; rather they are a few data points corresponding to nominal times when the target begins the next stage. For example, if the target has three stages, there will be three corresponding nominal staging times. Nominal staging times are used to determine the target's current stage number thisStage. The amount of elapsed time in the current stage tInstage is determined in a block 616 bytInstage=template_index−stageTimePast  (1)where stageTimePast is the nominal stage time of the previous stage, as illustrated in FIG. 7. In FIG. 7, the current time into the template is template_index. Block 618 of FIG. 6 represents the initialization of Kalman filter states and error covariances to produce the following equations (2), (3), (4), (5), and (6)X=XM  (2)is the initial position vector (set to measured value);{dot over (X)}={dot over (X)}M  (3)is the initial velocity vector (set to an initial estimate);K=1  (4)is the motor parameter (set to nominal value);delT=0  (5)is the initial estimate of error in template_index; anddelA=0  (6)is the initial estimate of error in angle of attack.
The TUBA Initialization function 514, represented as 600 of FIG. 6, includes a block 620 which represents the definition of the timing parameter template_index used to look up values from the nominal templateTemplate_index=(tInstage−delT)*K+stageTimePast  (7)From block 620, the logic of FIG. 6 flows to a block 622, which represents the definition of the time offset between the system time and template_indexoffsetError=template_index−tgtMeasTime  (8)where variable “tgtMeasTime” is the current system time, variously designated t or tm in FIG. 5.
TUBA uses a nine-state Extended Kalman filter which estimates the position and velocity vectors and three additional states. The three additional states are used to resolve the deficiencies associated with the use of a nominal acceleration profile. The TUBA filter equations are developed under the assumption that the target is either ballistic (falling under the force of gravity) or the specific force (such as thrust acceleration) is exactly known and can therefore be compensated for. It is also assumed that the target is not subject to significant atmospheric drag, which is reasonable in view of the high altitudes at which target tracking occurs. Alternatively, it is assumed that atmospheric drag can be properly compensated for. Equations (9), (10), and (11) model the target kinematics under these assumptions{umlaut over (X)}=Acc Gravity+Acc Centripetal+Acc Thrust+Acc Coriolis  (9)
                                          X            _                    ¨                =                                                            -                μ                            ⁢                                                          ⁢                              Z                _                                                                                                      Z                  _                                                            3                                -                                    ω              _                        ×                          (                                                ω                  _                                ×                                  Z                  _                                            )                                +                                                  Tacc                                      *                                          T                ^                            _                                -                      2            ⁢                          ω              _                        ×                                          X                _                            .                                                          (        10        )            Z=X+Re  (11)
where:
μ is the Earth gravitational constant, 3.986005*1014 m3/s2;
w is the magnitude of Earth's angular velocity, 7.29211574*10−5 rad/s;
|Tacc| is the boost acceleration magnitude; and
{circumflex over (T)} is the unit thrust vector.
The TUBA state vector is
                              s          _                =                  {                                                                      X                  _                                                                                                                          X                    .                                    _                                                                                    K                                                                    delT                                                                    delA                                              }                                    (        12        )            where:
X and {dot over (X)} are three-dimensional position and velocity vectors, respectively;
K is a scalar factor reflecting the target rocket motor, hot (K>1.0), cold (K<1.0), or nominal (K=1.0);
delT is the error in the initial time used to look up target parameters from the nominal templates; and
delA is the error in the estimate of the target's angle of attack.
The TUBA dynamics equations (i.e. the nonlinear TUBA state derivative equations) are
                                          s            _                    .                =                                            ⅆ                              ⅆ                t                                      ⁢                          s              _                                =                      {                                                                                                      X                      .                                        _                                                                                                                                                                                            -                          μ                                                ⁢                                                                                                  ⁢                                                  Z                          _                                                                                                                                                                            Z                            _                                                                                                    3                                                              -                                                                  ω                        _                                            ×                                              (                                                                              ω                            _                                                    ×                                                      Z                            _                                                                          )                                                              +                                                                                          Tacc                                                                    *                                                                        T                          _                                                ^                                                              -                                          2                      ⁢                                              ω                        _                                            ×                                                                        X                          _                                                .                                                                                                                                          0                                                                              0                                                                              0                                                      }                                              (        13        )            and are based on the assumed target kinematics set forth above. Additionally, it is assumed that K, delT, and delA are constants.
The thrust acceleration |Tacc|, is obtained from the nominal acceleration template using template_index and K estimate as|Tacc|=K*BoostAccLookup(template_index)  (14)
The unit thrust vector {circumflex over (T)} is calculated using the nominal angle of attack lookup, estimated position, velocity, and delA as
                    AoaEst        =                              AOALookup            ⁡                          (              template_index              )                                +          delA                                    (        15        )                                                      yL            ⁢                                                  ⁢            2                    _                =                              X                          .              ^                                _                                    (        16        )                                                      zL            ⁢                                                  ⁢            2                    _                =                              (                                                            X                                      .                    ^                                                  _                            ×                                                X                  ^                                _                                      )                    ×                                    yL              ⁢                                                          ⁢              2                        _                                              (        17        )                                                      T            _                    ^                =                                                            yL                ⁢                                                                  ⁢                2                            _                        *                          cos              ⁡                              (                AoaEst                )                                              -                                                    zL                ⁢                                                                  ⁢                2                            _                        *                          sin              ⁡                              (                AoaEst                )                                                                        (        18        )            where the template_index is defined by equation (7).
Referring once more to FIG. 5, the logic flows from initialization block 514 to a state propagation step 518, which is part of time update block 516. In the following description, a ^ notation is used to denote filter estimates of the respective variables, and a superscript (*)− represents a time update. Time propagation of the TUBA state vector in block 518 of FIG. 5 is performed by numerically integrating the state derivative vector from the previous time ti-1 to the current time ti where the subscript i refers to the filter cycle iteration
                                                                        s                _                            ^                        ⁡                          (                              t                i                            )                                -                =                                                            s                ^                            _                        ⁡                          (                              t                                  i                  -                  1                                            )                                +                                    ∫                              t                                  i                  -                  1                                                                              t                                      i                    -                    1                                                  +                                  Δ                  ⁢                                                                          ⁢                  t                                                      ⁢                                                                                                      s                      .                                        ^                                    _                                ⁡                                  (                  τ                  )                                            ⁢                                                          ⁢                              ⅆ                τ                                                                        (        19        )            A high order numerical integration algorithm such as the 2nd order Runge Kutta algorithm might be used for the integration process. The incremental time step, Δt, refers to either the nominal update cycle time or the incremental time step from the last cycle time to the current measurement time tM (i.e. Δt=tm−ti-1).
From state time propagation block 518, the logic of FIG. 5 flows to a block 520. Block 520 represents calculation of the state transition matrix for the extended TUBA Kalman filter algorithm. The Jacobian matrix corresponding to the state vector and dynamics described in equations (12) and (13) is
                    J        =                              [                                                  ⁢                                          ∂                                                      s                    .                                    _                                                            ∂                                  s                  _                                                      ]                    =                                    [                                                          ⁢                                                                                                                  ∂                                                                              x                            _                                                    .                                                                                            ∂                                                  s                          _                                                                                                                                                                                                        ∂                                                                              x                            _                                                    ¨                                                                                            ∂                                                  s                          _                                                                                                                                                                                                        ∂                                                  K                          .                                                                                            ∂                                                  s                          _                                                                                                                                                                                                                                  ∂                          del                                                ⁢                                                  T                          .                                                                                            ∂                                                  s                          _                                                                                                                                                                                                                                  ∂                          del                                                ⁢                                                  A                          .                                                                                            ∂                                                  s                          _                                                                                                                                ]                        =                          [                                                          ⁢                                                                                          0                                              3                        ×                        3                                                                                                                        I                                              3                        ×                        3                                                                                                  0                                                        0                                                        0                                                                                                                                      ∂                                                                              X                            _                                                    ¨                                                                                            ∂                                                  X                          _                                                                                                                                                                        ∂                                                                              X                            _                                                    ¨                                                                                            ∂                                                                              X                            _                                                    .                                                                                                                                                                        ∂                                                                              X                            _                                                    ¨                                                                                            ∂                        K                                                                                                                                                ∂                                                                              X                            _                                                    ¨                                                                                            ∂                        delT                                                                                                                                                ∂                                                                              X                            _                                                    ¨                                                                                            ∂                        delA                                                                                                                                                                                                                                                                                                                                                            0                                              1                        ×                        9                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                0                                              1                        ×                        9                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                0                                              1                        ×                        9                                                                                                                                                                                                                                                                                                            ⁢                                                          ]                                                          (        20        )            where components of the Jacobian corresponding to ∂{umlaut over (x)}/∂s are defined below in equations 21 to 26
                                          ∂                                          X                _                            ¨                                            ∂                          X              _                                      =                                                                              -                  μ                                ⁢                                                                  ⁢                                                      Z                    _                                    ^                                                                                                                                          Z                      ^                                        _                                                                    3                                      ⁡                          [                              I                                  3                  ×                  3                                            ]                                -                                    3                                                                                                            Z                      ^                                        _                                                                    2                                      ⁢                                                            Z                  ^                                _                            ·                                                                    Z                    ^                                    _                                ′                                              -                                    [                              [                                                      ω                    ^                                    _                                ]                            ]                        ·                          [                              [                                                      ω                    ^                                    _                                ]                            ]                                                          (        21        )                                                      ∂                                          X                _                            ¨                                            ∂                                          X                _                            .                                      =                              -            2                    ·                      [                          [                                                ω                  ^                                _                            ]                        ]                                              (        22        )            Note that the [[•]] notation denotes a skew symmetric matrix of the vector argument.
            ∂                        X          _                ¨                    ∂      K        ⁢          ⁢  and  ⁢          ⁢            ∂                        X          _                ¨                    ∂      delT      are performed numerically by choosing some small value for ∂K and ∂delT. The resulting equations are
                                          ∂                                          X                _                            ¨                                            ∂            K                          =                              a            ⁢                                                  ⁢            1                    +                                    K              ^                        *                                                            a                  ⁢                                                                          ⁢                  2                                -                                  a                  ⁢                                                                          ⁢                  1                                            dk                                                          (        23        )                                                      ∂                                          X                _                            ¨                                            ∂            delT                          =                              -                          K              ^                                *                                    (                                                a                  ⁢                                                                          ⁢                  3                                -                                  a                  ⁢                                                                          ⁢                  1                                            )                        dt                                              (        24        )            where:dk=0.001dt=0.01a1=BoostAccLookup(template_index)a2=BoostAccLookup(template_index+(tInstage+dt)*dk)a3=BoostAccLookup(template_index+dt)The partial derivative of acceleration with respect to the error in angle of attack,
            ∂                        X          _                ¨                    ∂      delA        ,is given by
                                          ∂                                          X                ¨                            _                                            ∂            delA                          =                              K            ^                    *                      BoostAccLookup            ⁡                          (              template_index              )                                *                      ∂                                          T                _                            ^                                                          (        25        )            where:∂{circumflex over (T)}=−yL2*sin(AoaEst)−zL2*cos(AoaEst)  (26)
The partial derivative of {circumflex over (T)} is taken with respect to angle of attack only, since variations of the thrust vector with respect to position and velocity are minimal. {circumflex over (T)}, AoaEst, yL2, and zL2 are defined in equations (18), (15), (16), and (17), respectively.
From Jacobian computation block 520 of FIG. 5 the logic flows to a block 522, representing the covariance propagation portion of the time update 516. The state transition matrix used for the time propagation of the TUBA error covariance can be approximated as
                    Φ        ≈                  I          +                      J            ⁢                                                  ⁢            Δ            ⁢                                                  ⁢            t                    +                                                    J                2                            ⁢              Δ              ⁢                                                          ⁢                              t                2                                      2                                              (        27        )            Time propagation of the TUBA error covariance matrix P is performed with the equationP(ti)−=ΦP(ti-1ΦT+Qi  (28)where:
Q is the 9×9 TUBA state noise matrix whose diagonal elements are chosen based on tuning considerations.
From covariance propagation portion or block 522 of time update 516 of FIG. 5, the logic flows to a Template Parameter Update block 524 and then to a Staging Logic Early Thrust Termination Block 542, both described in more detail below.
The logic of FIG. 5 arrives at a Gain Computation portion 528 of a Measurement Update block 526. The Kalman gain matrix is calculated in block 528 using the measurement matrix and the TUBA error covariance matrixK=P(ti)−·HT·(H·P(ti)−·HT+R)−1  (29)where:H=[I3×3 03×3 03×3]H is the measurement matrix, and R is the measurement noise covariance matrix associated with the currently reporting sensor.
The logic flows from block 528 of FIG. 5 to a State Measurement Update block 530. State Measurement Update block 530 computes the measurement residuals and updates the TUBA state vectorŝi=ŝi−+K·(Xm−H·{circumflex over (X)}−)  (30)
Finally, the logic of FIG. 5 reaches the Covariance Measurement Update function 532, which performs the measurement update of the state covariance matrixP(ti)=(I−K·H)·P(ti)−  (31)and the logic returns to time update block 516, with updated time T=Tm+Δt for the next calculation cycle, by way of a logic path 540.
Template Parameter Update block 524 of FIG. 5 updates the time index, template_index, used in the next filter update pass to look up the thrust acceleration and angle of attack. Template parameters tInstage and template_index are defined in equations (1) and (2), respectively, and were also described in conjunction with block 514 of FIG. 5. At each update after initialization, tInstage and template_index are updatedtInstage=tInstage+Δt  (32)where Δt is the measurement update intervaltemplate_index=(tInstage−del{circumflex over (T)})*{circumflex over (K)}+stageTimePast  (33)where stagetimepast is the nominal past stage.Note that the initial error in tInstage due to the initial guess at the value of target time after launch is removed in the calculation of template_index via the filter state del{circumflex over (T)}, and the acceleration profile variation due to motor differences is corrected by the filter state {circumflex over (K)}.
The TUBA target staging and early thrust termination logic block 542 of FIG. 5 provides timely estimates of the target's burn times and detects when early thrust termination has occurred. Block 542 is not mandatory for estimating the states of a boosting target using nominal template profiles of target (a) thrust acceleration, (b) altitude, (c) speed, and (d) angle-of-attack as a function of time after launch, updated using the states of a filter to correct for template indexing error, template boost acceleration variations, and template angle-of-attack variations. However, the Staging Logic Early Thrust Termination block 542 can improve the accuracy of the target state estimation by providing timely transition to the ballistic phase.
If the target is still thrusting, the estimated burn time for the current stage is calculated in block 542 of FIG. 5 according to equation (34)
  tBOEst  =                    (                  stageTimeCurrent          -          stageTimePast                )                    K        ^              +    stageTimePastEst    -    offsetError    +                  de        ⁢                  l          ^                ⁢        T                    K        ^            where:
stageTimeCurrent=nominal current stage burn out time
stageTimePast=nominal past stage burn out time
stageTimePastEst=estimated value of the previous stage burn out time
{circumflex over (K)} and de{circumflex over (l)}t=state estimates from the filter
offsetError=initial difference between template_index and measured time.
Staging times for subsequent stages are updated relative to the current stage time estimate.
To account for the uncertainties in the estimated burn out time, a staging window (e.g. 3 to 4 seconds) is set on either side of the estimated current stage burn out time in block 542 of FIG. 5. If the target is within the staging window, the estimated burn out time is used to determine if the target has entered the next stage or has entered the ballistic phase. When a staging event has occurred, filter estimates for {circumflex over (K)} and del{circumflex over (T)} are reset to their initial values. The target stage number is incremented, and template timing parameters are reinitialized as
{circumflex over (K)} = 1(35)del{circumflex over (T)} = 0thisStage = thisStage + 1offsetError = stageTimeEstPast − tgtMeasTime tInstage = tgtMeasTime − stageTimeEstPasttemplate_time = (tInstage − del{circumflex over (T)}) * {circumflex over (K)} + stageTimePastwhere:
tgtMeasTime is current system time.
When the target has exited a staging window, the filter error covariance for K and delA are reset to their initial default values by the Staging Logic portion of block 542 of FIG. 5. When the target is ballistic, the error covariance for K, delT, and delA are all reinitialized to some small value, (e.g. 0.00000001).
Liquid-propellant rocket engines control the thrust by varying the amount of propellant that enters the combustion chamber. By stopping the flow of propellants into the combustion chamber, these targets can effectively terminate their thrust prior to the nominal end of boost time. Solid-propellant rockets are not as easy to control as liquid rockets. Once started, the propellants burn until they are gone. However, some solid-fuel engines have hatches on their sides that can be cut loose by remote control to release the chamber pressure and terminate thrust. Early thrust termination poses a particularly difficult problem for tracking since it is an unknown change in acceleration applied at an unknown time. If the filter continues to assume nominal boost acceleration when in actuality the rocket has early thrust terminated, potentially huge errors would result in the estimated states which might well render invalid any fire control solution. The early thrust termination logic portion of block 542 of FIG. 5 attempts to overcome this difficulty. More particularly, at 20 seconds (or at some other selected time) prior to the estimated final burn out time, TUBA increases the K covariance values and observes for changes in K estimates. At this point, the K factor has settled to the correct value reflecting the state of the rocket motor. Therefore, a sudden and consistent change in K can only be caused by a substantial change in the target's acceleration. Such a sudden, consistent change in K indicates that the target has early thrust terminated, and gives an indication of the time of thrust termination which can be used in the overall state estimates.
Thus, a method according to the Luu et al. patent application is for estimating the states of a boosting target using nominal templates (300). The templates (300) comprise nominal profiles of the target (a) thrust acceleration, (b) altitude, (c) speed, and (d) angle-of-attack as a function of time after launch, all updated using the states of a filter to correct for template indexing error, template boost acceleration variations, and template angle-of-attack variations, thereby producing estimated target states unencumbered by motor burn boost variations and providing for accurate determination of end-of-boost.
According to another aspect of the method of the Luu et al. patent application, for estimating at least the position and velocity state of a boost vehicle target, the state estimates include the position and velocity, and at least one of (a) variation in engine burn rate, (b) error in time after launch and (c) error in angle of attack. The method comprises the step of sensing at least the position of the target to thereby generate target-position information. The method also comprises the step of providing a template (300) characterizing nominal values of acceleration of the boost vehicle target due to engine thrust, speed, altitude, and angle-of-attack, all as a function of time after target launch. The target altitude is calculated or determined, possibly from the target-position information, to generate calculated target altitude. A time index (Io) is established by entering the template (300) at the calculated target altitude (Alto). At least that (or those) nominal values of time after launch (TAL), acceleration (ABx), speed (Vx), and angle of attack (AOAx) corresponding to the calculated target altitude are read from the template, where the subscript x represents any value. Both the state estimate and the associated covariances of an extended Kalman filter are initialized with the nominal value of time after launch, acceleration, speed, and angle of attack corresponding to the calculated target altitude. The state estimates and covariances of the Kalman filter are propagated in time (516) to generate a time-updated vector of state estimates representing the predicted state of the target, using the template profile information (524) updated by at least one of (a) variation in engine burn rate, (b) error in time after launch and (c) error in angle of attack. The states are updated (526) using the measurements, to thereby produce a vector of measurement updated state estimates and covariances of the boost vehicle target, including position and velocity, and at least one of (a) variation in engine burn rate, (b) error in time after launch and (c) error in angle of attack.
The TUBA method for determining the state of a target relies, at least in part, on identification of the target, which is to say identification of the type of missile (and therefore the parameters of the missile). This reliance also applies to other template-based missile tracking methods. The identification of the target relies, at least in part, on the matching of the sensed parameters to the nominal templates after the filters have settled. However, the sensed parameters of a given missile may deviate from the nominal, and possibly overlap with the nominal values of other missile types. Thus, it is conceivable that the wrong nominal filter will be selected as the best representation of the target and its parameters. This, in turn, can be expected to result in less-than-optimal predictions being used in an antimissile fire-control system. Improved or alternative methods for identification of the missile type are desired to aid in establishing its parameters for tracking.